Prom waves to avalanches: two different mechanisms of sandpile dynamics 
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Time series resulting from wave decomposition show the 
existence of different correlation patterns for avalanche dy- 
namics. For the d — 2 Bak-Tang-Wiesenfeld model, long 
range correlations determine a modification of the wave size 
distribution under coarse graining in time, and multifractal 
scaling for avalanches. In the Manna model, the distribution 
of avalanches coincides with that of waves, which are uncor- 
rected and obey finite size scaling, a result expected also for 
the d = 3 Bak et al. model. 
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Full information on a self-organized critical (SOC) pro- 
cess is contained in the time series, if the time step is 
the most microscopic conceivable. The self-similarity of 
such a process, due to its intermittent, avalanche charac- 
ter, should be revealed by the scaling of the time auto- 
correlation or its power spectrum. In spite of this, time 
series analyses have been seldom performed on sandpiles 
or similar systems, and mostly concerned the response 
to finite, random external disturbances, i.e. the prob- 
lem of 1/f noise [Q. Most efforts in the characterization 
of SOC scaling concentrated on probability distribution 
functions (PDF's) of global properties of the avalanche 
events, which occupy large intervals of the microscopic 
evolution time ||. The numerical analysis of such PDF's 
is often difficult, and universality issues can not be easily 
solved. The situation is even more problematic if, like for 
the two-dimensional (2D) Bak-Tang-Wiesenfeld sandpile 
(BTW) 0], the usually assumed finite size scaling (FSS) 
form reveals inadequate for the PDF's, and needs to be 
replaced by a multifractal one pj . These findings rise the 
additional issue of why the 2D BTW displays such mul- 
tifractal scaling, while apparently very similar sandpiles, 
like the Manna model (M) §], do not HH]. 

In recent theoretical approaches to the BTW and other 
Abelian sandpiles (8), a major role has been played by 
the waves of toppling into which avalanches can be de- 
composed ||. For the BTW, the PDF's of waves, as 
sampled from a large collection of successive avalanches, 
obey FSS with exactly known exponents |l0|. By ana- 
lyzing the succession of waves as a stationary time series, 
one could hope to determine the statistical properties of 
avalanches. However, so far, no precise information on 
the correlations of such series was obtained. 

In this Letter we generalize the wave description to 
the M in 2D. The study of the respective wave time 



series reveals that BTW and M are prototypes of two 
very different scenarios for avalanche dynamics. In the 
M case, successive wave sizes are totally uncorrelated. As 
a consequence, avalanche and wave PDF's have identical 
scaling properties, consistent with FSS. For the BTW, to 
the contrary, wave sizes show long range correlations and 
persistency in time. PDF's of "block variables" given by 
sums of n successive wave sizes show that these correla- 
tions are responsible for the fact that avalanche scaling 
differs from the wave one and has multifractal features. 
For the 3D BTW, on the other hand, our results suggest 
validity of a scenario identical to the M one, and lead 
to conjecture exact avalanche exponents coinciding with 
those of the wave PDF. 

The 2D BTW § is defined on a square L x L lattice. 
An integer Zi, the number of grains, is assigned to site 
i. Starting from a stable configuration (zi < z c = 4, 
Vi), a grain is added at a randomly chosen site; after 
each addition, all sites exceeding the stability threshold, 
Zk > z c , undergo toppling, distributing one grain to each 
one of the nearest neighbors. The topplings, which dis- 
sipate grains when occurring at the edges, continue un- 
til all sites are stable, and a new grain is added. The 
s topplings between two consecutive additions form an 
avalanche. After many additions, the system organizes 
into a stationary critical state. 

Manna || studied a two-state version, M, of the sand- 
pile. The sites can be either empty or occupied; grains are 
added randomly, and when one of them drops onto an oc- 
cupied site, a "hard core" repulsion pushes two particles 
out to randomly chosen nearest neighbors. Compared to 
the BTW, in which toppling is deterministic, this model 
has an extra stochastic ingredient in the microscopic evo- 
lution, and z c — 1. 

One can define a special ordering in the topplings 
of the BTW by introducing the wave decomposition of 
avalanches 0]. After the site of addition, O, topples, any 
other unstable site is allowed to topple except possibly O. 
This is the first wave of toppling. If O is still unstable, 
it is allowed to topple once again, leading to the sec- 
ond wave. This continues until O becomes stable. Thus, 
an avalanche is broken into a sequence of waves. Dur- 
ing a wave, sites topple only once, and for an m-wave 
avalanche s = X)feLi s fe; where s& are the topplings of 
the fc-th wave. Following the definitions for the BTW, 
we implement here a wave decomposition of avalanches 
also for the M. Unlike for the BTW, sites involved in 
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a wave may topple more than once. Furthermore, the 
toppling order chosen implies now a peculiar sequence of 
stable configurations visited upon addition of grains, but 
the realization probabilities of possible configurations are 
independent of this order || . 

The BTW wave size PDF has FSS form P w (s) ~ 
s~ T ™p w {s/ L Dm ), with t w = 1 , D w = 2 in 2D, and p w 
a suitable scaling function HftO). For the M waves we 
also found a PDF obeying FSS, with r w = 1.31 ± 0.02 
and D w = 2.75 ± 0.01. These exponents are remark- 
ably consistent with those estimated for avalanches []5||| . 
Such a coincidence certainly does not apply to the BTW. 
Attempts to derive exact BTW avalanche exponents 
within FSS were based on the observation that often 
waves within an avalanche show a final, long contrac- 
tion phase, and on a scaling assumption for the corre- 
sponding Sk — s/c+i JTlj - A more adequate ansatz for 
the conditional PDF of Sk+i given Sk |l2|], and Marko- 
vianity assumptions, did not help in better characterizing 
avalanche scaling along these lines |0|. In fact the 2D 
BTW obeys a multifractal form of scaling, which is not 
catched by such simplified schemes 0j. 

The wave time series {sfc} provide coarse grained dy- 
namical descriptions. In the L — > oo limit, these de- 
scriptions are infinitely rescaled with respect to those at 
microscopic time scale, but still infinitely finer than the 
mere records of successive avalanche sizes. This interme- 
diate time scale reveals essential in order to understand 
the dynamics indside avalanches, whose size sequence we 
found to be uncorrelated in the sense discussed below for 
waves. The microscopic scale gave no significant results 
in comparing the two models, since, at that level, simi- 
lar strong correlations exists in both of them, due to the 
parallel updating algorithm. 

We determined for BTW and M the autocorrelation 



C(t,L) = 



(s k+t s k ) L - [l 2 

(4) L 



(1) 



with t — 1,2,..., and /i = (sk) L , the time averages being 
taken over up to 10 7 waves for L = 128, 256, 512, 1024 
and 2048. A first, striking result is that, as soon as t > 
(C(0, L) — 1 by normalization), waves are uncorrelated 
in the M. Indeed, as L grows, C manifestly approaches 
as soon as t > (Fig. |l|). To the contrary, C is long 
range for the BTW, because it approaches only for t 
exceeding the maximum number of waves in an avalanche 
(Fig. Q), which we found to scale ~ L, for L — > oo. For 
the BTW we further tested the FSS form 

C(t,L)=t- T 'g(t/L D ") (2) 

on the basis of the L — > oo scaling of the moments ijllpl 



(?) 



(3) 



Cc(<?) 



D c (q-T c + l) if?>r c -l, 
if q < t c - 1. 



(4) 



Fig. H shows the extrapolated cr c (q), which has an ap- 
proximately linear part for 1 < q < 4, consistent with 
D c = 1.02 ± 0.05 and t c = 0.40 ± 0.05. The curvature 
for q < 1 is due to the fact that, for finite L, a loga- 
rithm (t c = 0) can not be easily distinguished from a 
power law with r c > 0. The inset of Fig. [2] shows this 
logarithmic growth of (t q ) L . We conclude that, for the 
BTW, a simple power law tail C(t, oo) ~ t~ Ta , is a first, 
rough, approximation. The increment y(t) = Sk 
is comparable to the trail of a fractional Brownian mo- 
tion with Hurst exponent H = (2 — t c )/2 [flBI , such that 
H = 0.80 ± 0.03 for the BTW. H can be measured di- 
rectly from the fluctuation 



F(t,L) = (Ay(tf) L - (Ay(i)) J 



1/2 



(5) 



y(k 

4-H 



t) — y{k) p5|, which should scale 



with Ay(t) 

as F(t, oo) ~ t H . Fig. || reports tests of this scaling for 
L = 2048. For the BTW, H - 0.85 at low t, in agreement 
with t c — 0.40±0.05. A crossover to H = 1/2 is observed 
for large t. The crossover time of course increases with L. 
H = 1/2 corresponds to a process with C exponentially 
decaying with t, or with C(i, oo) = S t .o, as is the case for 
the M. Thus, the crossover is due to the fact that, beyond 
the maximal time duration of avalanches, waves are un- 
correlated. H ^ 1/2 implies long range correlations, as 
we find for the BTW. Furthermore, H > 1/2 corresponds 
to persistency: an increasing or decreasing trend of y in 
the past mostly implies a similar tendency in the future. 
This accounts for the observed expansion and contrac- 
tion phases in avalanche growth |ll|] . The above features 
of C should be responsible for the peculiar long range 
on-site correlations of the noise expected when mapping 
the BTW onto a discrete interface growth equation [jl6| . 

The BTW has been shown to display a non constant 
gap in the high q moments of some avalanche PDF's 0. 
Thus, it should not surprise if, unlike assumed in Eq. (|), 
also C would manifest similar multifractal properties. In- 
deed, a more accurate analysis shows that, for high g, the 
gap da c (q) jdq grows slowly with q, beyond the above esti- 
mate of D c . This multifractal character can be embodied 
in the more general scaling ansatz M 



C(t = L a ,L)=L^ a > 



(L » 1) 



(6) 



FSS would imply the piecewise linear form 



with a nonlinear singularity spectrum f c > — oo in an 
a-interval covering the whole range of possible gaps and 
linked to a c by Legendre transform. C would satisfy 
the FFS ansatz (||) only if f c were a linear function, i.e. 
f c (a) = -(r c - l)a if a G [0,D C ], f c (a) = -oo other- 
wise. The FSS picture given above is in fact only an 
approximation. This explains also the slight curvature 
of the F plot for low t, which makes the direct measure- 
ment of H ambiguous (Fig. ^|). Rather than attempting 
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a more precise determination of a c and f c , below we clar- 
ify the difference between BTW and M, and the origin 
of multiscaling in the former, in the light of probabilistic 
concepts. 

Waves have a relatively simple behavior. So, in a renor- 
malization group (RG) spirit, we can coarse grain the 
time, by looking at the PDF, p( n )(s, L), of the sum of the 
sizes of n consecutive waves, regardless of the avalanche 
they belong to. Since avalanches are constituted by an 
infinite number [ [l7f of waves for L — > oo, by sending 
also n — > oo, we expect p( n > to approach the PDF of 
avalanche sizes as a RG fixed point. This approach can 
be monitored on the basis of the effective moment scaling 
exponents defined by 

<*'} n , z = J dss«p( n \s,L)~L°^. (7) 

For the M, <7 Si „(g) does not depend on n, and, within 
our accuracy, is equal from the start to cr s (q), such that 
/ ds s q P(s,L) ~ L as{q \ with P(s,L) representing the 
avalanche size PDF. To the contrary, in the BTW, as n 
increases, a s ^ n {q) varies and moves gradually towards the 
appropriate a s (q) (Fig. |]). 

The result for the M can be explained on the basis of 
the fact that PDF's of independent variables, satisfying 
a scaling ansatz of the type (^|), have a spectrum which 
does not change under convolution pq |. For example, 
if we put P w (s = L a , L) ~ £/«.(«)-« w here f w is the 
spectrum of the wave size PDF, one can verify that 

pW(L a ,L)= (8) 

J d(lf) P w (L a - LP,L)P W (LP,L) ~ . 

Thus, also the spectrum associated with pW(s,L), 
which is the convolution of n P w 's, does not depend on 
n. This implies the n-independence of a s ^ n (q), which is 
determined once f w is given. Now, let P(s, to, L) be the 
probability of having an avalanche with s topplings and 
to waves: due to the uncorrelated character of different 
wave sizes, for the M P(s,m,L) is also the convolution 
of to P w 's, i.e P{s,m,L) = P (m ^(s,L). Therefore, if 
P(m, L) is the PDF of the total number, to, of waves in 
an avalanche, one has 

P(s,L)=J2P^(s,L)P(m,L), (9) 
m 

such that 

L CT ^) ~ J2 L a °-™ (q) P{m, L), (fO) 

m 

which, together with the above results, imply <r Sin (q) = 
<y s {q)- The coarse graining of waves in the M does 
not modify the block PDF spectrum, which is from the 
start at its fixed point, representing also the scaling of 
avalanches. 



In the BTW case, the nontrivial composition of cor- 
related waves is responsible for the change with n of 
the effective singularity spectrum of P^ (Fig. ^). For 
n = 1, the simple linear FSS form f w (a) = —(t w — l)a 
(a G [0,-D,,;]) applies, while for n — > oo one should re- 
cover the nonlinear form needed to describe the com- 
plex scaling of P(s,L) Q. In practice, sampling lim- 
itations prevent us from reaching very high n. How- 
ever, even if convergence is relatively slow, the ten- 
dency of (J s ,n{q) to move towards u s (q) is very manifest. 
We could clearly detect the increase with n of the gap 
d<J s ,n{q) / dq, at fixed high q: for example, we determined 
da s [ n (4)/dq ~ 2.00,2.20,2.25,2.38, for n = 1,8,12,24, 
respectively. Furthermore, while the asymtoptic gap 
for P w (s,L) gets readily to the maximum value, i.e. 
dcr s .i(q) I 'dq = 2 for q > 1, as soon as n > la con- 
stant gap could not be detected for <7 Sj71 (g); this confirms 
the progressive appearance of effective multifractal scal- 
ing for the block variable. 

In summary, a comparative analysis of wave time se- 
ries shows that the different forms of universal scaling 
in the BTW and M are determined by distinct dynam- 
ical correlation patterns. For the M case, the wave 
level of description is in fact coarse grained enough to 
fully account, without further modifications, also for the 
avalanche level. Waves are uncorrelated and their PDF, 
satisfying FSS, coincides with the avalanche PDF, as far 
as exponents are concerned. To the contrary, in the 
BTW, under coarse graining, long time correlations sub- 
stantially modify the scaling properties of waves, deter- 
mining also multiscaling features. We regard these as 
prototype mechanisms to be generally expected in sand- 
pile and similar SOC models. The 2D BTW behavior is 
probably less generic than the M one. For example, in the 
3D BTW, the fact that avalanches are constituted by an 
average number of waves which remains finite for L — > oo, 
strongly suggests a M type mechanism and FSS for the 
avalanche PDF. Indeed, even if subsequent waves are not 
strictly uncorrelated, like in the 2D Manna model, C 
should decay extremely fast, with an L independent time 
cutoff. Thus, coarse graining in time can not substan- 
tially modify the PDF of block variables, with respect to 
P w . We expect the exactly known wave exponents to ap- 
ply also to the avalanche size PDF. Indeed, numerically 
determined avalanche exponents for 3D BTW turn out 
to be strikingly close to the wave values pp| . 
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FIG. 1. For the BTW model, as L grows, C decays more 
and more slowly, (L = 128, 256, 512, 1024, 2048 from left to 
right). C(t > 1, L) ~ for the M (inset). 
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FIG. 2. The dashed curve interpolates o" c (g) for 1 < q < 4; 
it has slope D = 1.02 and intercept 0.62. The inset shows 
semi-logarithmic plots of (t q ) L vs log 2 L for q = —1, —0.5. 
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FIG. 3. Plots of F(t,L = 2048) for the BTW (continuos) 
and the M (dashed, slope = 0.5). 
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FIG. 4. a a>n {q) for the BTW, n = 1,8,24 (open circles, 
from the bottom). The full circles show er s (g). 



